arXiv:1505.02959vl [gr-qc] 12 May 2015 


Expanding universe with nonlinear 
gravitational waves 

Taishi Ikeda*t Chul-Moon Yoo^^, and Yasusada Nambu^^ 

^Departure of Physics, Graduate School of Science, Nagoya 
University, Nagoya 464-6602, Japan 

Abstract 

We test the validity of Isaacson’s formula which states that high 
frequency and low amplitude gravitational waves behave as a radia¬ 
tion fluid on average. For this purpose, we numerically construct a 
solution of the vacuum Einstein equations which contains nonlinear 
standing gravitational waves. The solution is constructed in a cubic 
box with periodic boundary conditions. The time evolution is solved 
in a gauge in which the trace of the extrinsic curvature K of the time 
slice becomes spatially uniform. Then, the Hubble expansion rate H 
is defined by H = — A/S and compared with the effective scale factor 
L defined by the proper volume, area and length of the cubic box. 

We find that, even when the wave length of the gravitational waves 
is comparable to the Hubble scale, the deviation from Isaacson’s for¬ 
mula H oc L~^ is at most 3% without taking a temporal average and 
is below 0.1% with a temporal average. 


1 Introduction 

The global aspect of our universe is well described by a homogeneous and 
isotropic FLRW(Friedmann-Lemaitre-Robertson-Walker) metric. The homo¬ 
geneity over the scale of lOOMpc is well established by many observations. 
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At the same time, local inhomogeneity of our universe provides us a lot of 
information. One of the classical issue involved with the inhomogeneity is 
the so called backreaction problem(see e.g. [l]-0]). Evaluation of the backre- 
action due to the local inhomogeneity on the global expansion law has been 
attracted much attention because large backreaction may cause significant 
systematic error for evaluation of the energy components of our universe. 

One typical example of treating the backreaction is Isaacson’s formula 0 
E]. For gravitational waves with high frequency and low amplitude, Isaacson’s 
formula provides the coarse-grained effective energy momentum tensor which 
satishes the traceless condition. Since the energy momentum tensor which is 
compatible with the homogeneity and isotropy is given by the perfect fluid 
form, Isaacson’s formula states that the high frequency and low amplitude 
gravitational waves behave like a radiation fluid on average (see Appendix 
[0. Therefore, for a spatially flat FLRW universe with short wavelength 
gravitational waves, the relation between the Hubble parameter H and the 
scale factor a is expected to be if oc a~^. 

The aim of this paper is to test the validity of Isaacson’s formula beyond 
the high frequency and low amplitude approximation. For the purpose, we 
consider the inhomogeneous universe which only contains the gravitational 
waves and calculate the time evolution of this universe. Because this universe 
is assumed to have no symmetry, it is inevitable to rely on a numerical 
computation to obtain its dynamics and we apply numerical relativity to 
tackle this problem. Over the past 20 years, numerical relativity has been 
extensively developed mainly for isolated systems. Recently, several authors 
applied the numerical relativity to the expanding or contracting universe 
models with aligned black holes and discussed the backreaction problem [7^ 

[in]. 

In this paper, we investigate the validity of Isaacson’s formula by solving 
the full Einstein equations with the use of numerical relativity. First, we 
construct the initial data which contains nonlinear gravitational waves by 
solving constraint equations of vacuum Einstein equations. This initial data 
is prepared in a cubic box with periodic boundary conditions between each 
pair of opposite faces of the box and corresponds to standing wave modes in 
the linear approximation. Concerning the time evolution, the backreaction 
effect of gravitational waves causes expansion or contraction of this universe. 
We use the gauge condition in which the trace of the extrinsic curvature K is 
spatially uniform. Then, the Hubble parameter is defined by if = —K/3. We 
introduce the effective scale factors L from the proper volume, proper area. 
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and proper length of this box. Isaacson’s formula states that time evolution 
of this system is same as the universe with the radiation fluid in the short 
wavelength region and the Hubble parameter and the effective scale factors 
obey the relation 



( 1 ) 


where 6 is a constant. We evaluate the relation between the Hubble parameter 
H and the effective scale factor L in the numerically generated inhomoge¬ 
neous universe and test the validity of this formula. 

This paper is organized as follows. In Sec. 2, the initial data are con¬ 
structed by solving constraint equations of the vacuum Einstein equations. 
We describe the time evolution of our model in Sec. 3 and the relation be¬ 
tween the Hubble parameter H and the effective scale factor L is investigated. 
Furthermore, we report dependence of the initial amplitude on the dynamics 
in superhorizon scale and the effect of temporal average of the volume. Sec. 4 
is devoted to a summary and discussion. We use units in which c = G = 1 
throughout this paper. 


2 Set up of initial data 

2.1 Construction of initial data 

As is mentioned in the introduction, we numerically construct a solution of 
the vacuum Einstein equations, which contains nonlinear standing gravita¬ 
tional waves in the cubic box with periodic boundary conditions. Naively 
thinking, since the gravitational waves have a positive energy, it causes the 
gravitational attractive force, so this system is unstable. Actually, as we will 
see later, this universe inevitably expands or contracts. In this section, we 
describe how to construct the initial data. 

The initial data set is generated by solving the following Hamiltonian and 
momentum constraints: 


n + K"^- = 0 , ( 2 ) 

DjiG, - D,K = 0, (3) 

where TZ is the scalar curvature of 3-metric 'jij, Kij is the extrinsic curvature, 
K = 7 jjA'*A and Di is the covariant derivative associated with 7 ^^. For con¬ 
venience, 74 j and Kij are decomposed into the conformal factor T, conformal 
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3-metric 7 *^, and conformal traceless part of the extrinsic curvature 
follows: 


7 y = det(%) = 1 , 

Kij = (^Aij + l^ijK^ , = 0 . 


Aij as 


(4) 


Using these variables, the Hamiltonian constraint and the momentum con¬ 
straints are written as 


D.m - = 0 , (5) 

DjA^^ + QA^Wj In $ - -&K = 0, ( 6 ) 

3 

where Di and R are the covariant derivative and the scalar curvature asso¬ 
ciated with 7 jj, respectively. In order to construct the initial data, we solve 
these constraints with appropriate ansatzes. 

We assume that K is spatially constant and Aij = 0. These assump¬ 
tions make momentum constraints trivial. Furthermore, the Hamiltonian 
constraint is reduced to 

DiD^^ - = 0. (7) 

We also assume that the conformal metric has the following form: 

, 7 + ftP) i + fcw i + ;iP)\ 

( 8 ) 

where = j\{^) cos(a;oa:’^4 -g = x, y, z) [ 11 ]. A constant ojq specifies 

the coordinate wave number of the gravitational waves and 0^4 jg the phase 
of the gravitational waves. As we will see later, in the linear approximation, 
these ansatzes lead to a solution for standing gravitational waves with the 
amplitude ^^4 in tpjg paper, for simplicity, we assume 

= j^{v) = ( 9 ) 


By this assumption, three spatial axes are equivalent to each other. 

We consider periodic boundary conditions for each pair of opposite faces 
of the numerical box. The coordinate length of the edge of the box is set to 
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be Ao = 2ti/ojq. 0 Under this boundary condition, without loss of generality, 
we can set 

_ 0(y) _ 0(d = 0. (10) 

By integrating the Hamiltonian constraint over the box, we obtain 


K 


± 




( 11 ) 


where the hrst term in Eq. ([^ has been reduced to the boundary integral 
and eliminated due to the periodic boundary conditions. Since the trace part 
of the extrinsic curvature gives the volume contraction rate, our solution 
describes expanding or contracting universe corresponding to negative or 
positive value of the extrinsic curvature, respectively. 

Hamiltonian constraint ([7]) is invariant under the following scaling with 
a constant C\ 

^ , K ^ KIC‘^. (12) 

This ambiguity corresponds to a choice of the unit of the scale. In this paper, 
the normalization is hxed by 


/ = 1. (13) 

J box 

With the conditions fITT]) and (ir3|) . Eq. ([7]) is numerically integrated by using 
the successive over-relaxation (SOR) method. 

We show that our initial data contains gravitational waves in the linear 
approximation Al <C 1. It is worthwhile to note that, since R = 0{A^), we 
obtain K = 0{A) from Eq. (ITTD and T = 1 -|- 0{A?) from the Hamiltonian 
constraint (j7]). In this approximation, the conformal metric becomes 

~ 6ij + diag(h*^^^ — (14) 

and the fluctuation part jij — 5ij satishes the transverse traceless condition. 
Furthermore, for = 0, as the equation of motion of is ^ 7 *^ = — 

Aij = 0 means ^ 7 p = 0. Hence, the linearized form of our initial data 7 jj 
corresponds to the standing gravitational waves at the moment of maximum 
amplitude. 

Hn practice, 1 /8 region of the box with reflection boundary condition is enough because 
of the discrete symmetry. 


5 





Remaning parameter A corresponds to the initial amplitude of gravita¬ 
tional waves in the linear regime. In the short wavelength, the linear grav¬ 
itational waves have only oscillatory modes and their amplitude decreases 
(increases) with expansion (contraction) of the universe. On the other hand, 
for the super horizon wavelength, the two modes of the gravitational waves 
correspond to the decaying mode and the growing mode. The ratio between 
these two modes is fixed by the phase of the gravitational waves at the horizon 
crossing. In our initial data, the phase of the gravitational standing waves 
is fixed so that the waves have the maximum amplitude at the initial time 
and the change of the value A causes the change of the initial Hubble scale. 
Therefore the phase at the horizon crossing depends on A and the behavior 
of the gravitational waves in the long wavelength regime also depends on the 
initial amplitude A. We use initial data with A = 0.07,0.08,0.09,0.1,0.11 
and analyze the behavior of the metric both in the short wavelength region 
and long wavelength region. 

3 Time evolution 

To simulate the time evolution, we use the COSMOS code by appropriately 
modifying it for our purpose. The COSMOS code is an Einstein equation 
solver written in C-I--I- by means of BSSN formalism [T^ITB] . The algorithm 
of this code is similar to the SACRA code [H]. The COSMOS code has been 
used in papers [TOllTS] . In the original code, the 4-th order hnite differencing 
in space with uni-grid and the 4th order time integration with Runge-Kutta 
method in Cartesian coordinates are adopted. In this paper, we adopt the 
2nd order central differencing in space and the leapfrog method with a time 
filtering for time integration. Because of the time reversal of the Einstein 
equation, the time evolution is simulated not only for the expanding universe 
with K < 0, but also for the contracting direction from the same initial data 
with K > 0. Our simulation has been terminated when the violation of the 
Hamiltonian constraint exceeds 1%. 

3.1 Gauge condition 

Before performing the time integration, we need to fix the gauge conditions 
for the lapse function a and the shift vector /9®. In regard to the shift vector, 
we set /3* = 0 for simplicity. As for the lapse function, we use the “uniform 
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K gauge” condition which keeps spatially uniform K on each time slice. 

Let us derive the equation for the lapse function with the uniform K 
gauge. The time evolution equation of K is given by 

= -D,D^a + a . (15) 

The uniform K gauge condition implies that the left hand side in this equa¬ 
tion is spatially constant. Integrating this equation over the box, the time 
derivative of K is obtained by 

dK + , 

where the hrst term in the right hand side of Eq. (IT^ vanishes by virtue of 
Gauss’s theorem and the periodic boundary conditions. At every time step, 
we solve Eq. flTbll by the SOR method to determined a. 


3.2 Physical quantities 

In order to test Isaacson’s formula o. we investigate the relation between 
the effective Hubble parameter and the effective scale factor. Since our model 
is inhomogeneous, there is no unique dehnition of the Hubble parameter and 
the scale factor. Nevertheless, by virtue of the uniform K gauge condition, 
the Hubble parameter can be naturally dehned as 


H = - 


K 


(17) 


This dehnition coincides with the standard dehnition of H for the homoge¬ 
neous and isotropic universe model. 

One of the simplest dehnition of the scale factor is dehned from the proper 
volume of the box: 


Tv = 


dxdydz^ 


' box 


1/3 


(18) 


We can also dehne other ehective scale factors from the edge proper length 
and the surface proper area. These dehnitions are used in the several previous 
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Figure 1: The box with periodic boundaries. The definitions of the lines 0-8 and surfaces 
80 ( 2 : = 0, blue surface), Si(z = Ao/4, red surface) and 82 ( 2 : = Ao/2, green surface) are 
shown. 


works [TlISl ITOllTSHTT] . The effective scale factors from the proper length and 
the proper area are 


I^h(x,y) 


La{z) 



(19) 

( 20 ) 


Since Lp and La depend on x, y, z, we pick up the several characteristic 
positions as shown in Fig.[Tl La{z) is evaluated at z = 0, Ao/4, Ao/2 which are 
labeled as So, Si and S 2 , respectively. Li^{x,y) is evaluated at {x,y) = (0,0), 
(0, Ao/4),(0, Ao/2),(Ao/4,0) ,(Ao/4, Ao/4),(Ao/4, Ao/2),(Ao/2,0) ,(Ao/2, Ao/4), 
(Ao/2, Ao/2), which are labeled as 0-8 lines, respectively. 

In order to compare our numerical results with Isaacson’s formula ([T]) 
in a precise sense, we must take not only spatial average but also temporal 
average. We consider the spatial average is taken by the dehnition of effec¬ 
tive scale factors flT8|) - fl20|) . However, for temporal average, we do not have 
any appropriate scale for the average, since our model has no exact period¬ 
icity along the time direction. One possible way for the temporal average is 
proposed and discussed in the Sec J4.3[ 

We evaluate the difference between the effective scale factor from our 
numerical result and that of Isaacson’s formula as a function of the Hubble 
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parameter. The deviation is evaluated by 




T T Isa 

T Isa 

-^v 


( 21 ) 


T T Isa 

V'" - 

-^A 

SL{H,x,y)= . (23) 

where Ly^^H), Vl^{H,z) and L^^{H,x,y) represent Isaacson’s formula ([I]) 
and the coefficient b is determined by the least squares fitting of our numerical 
result in the region Ly < Aq, which is attained by time evolution along the 
expanding temporal direction. 


3.3 Convergence check 

Let us consider the value of Sy as a function of the grid spacing A. Since 
our numerical code has the second order accuracy, (5y(A) and its true value 
(5vtrue = <5v(0) are supposed to have the relation 

5v(A) = <5vtrue + dA^ + 0{A^). (24) 


where a coefficient d is determined by numerical results. Using the two data 
sets taken with A = Ai and A 2 , (5vtrue and d can be obtained as 


d 


JVtrue 


(^v(Ai) - 6y{A2) 

Af - A| ^ ' 

6y{A2)Al - 5 y{Ai)Al 

A? - A2 ^ 


(25) 

(26) 


For Ai = Ao/120 and A 2 = Aq/IOO, we calculate the value of d and ^vtme- 
The value of convergence (5v(A) — 5vtrue)/d is compared with A^ for each 
value of A in Fig. |2j We can clearly confirm that the second order convergence 
is achieved from Fig. [2l In the following, we use the numerical result with 
A = Ao/120. 
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Figure 2: The convergence of (i5v(A) — 6 vtrue)/d at selected time steps, ^vtrue and d are 
fixed by Eqs. (|25|) and (|26l) . 

4 Results of simulation 

4.1 Relation between Hubble parameter and effective 
scale factor 

In this subsection, we explain the relation between the Hubble parameter 
and the effective scale factors Ly, La, Li^ for the initial amplitude A = 0.1. 

Before discussing the result, we introduce the characteristic Hubble pa¬ 
rameter Hi o, i/o .5 and Hq i by 

= H-\ Li-(ilo.5) = 0-5H-\ Li-(hro.i) = (27) 

where = ^Jb/H. That is, corresponds to the horizon 

crossing Hubble time. The parameter b in each is determined by 

htting the numerical data obtained by time evolution in expanding temporal 
direction. We dehne the short wavelength region and the long 

wavelength region > H~^ in comparison with the Hubble scale. 

The behavior of Ly{H) and Sy{H) is shown in Fig. [3l 
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Figure 3: The behavior of Ly{H). Isaacson’s formula is Ly^^H) = y/b^/H with by = 
0.224. 


These figures show the oscillation of the scale factor which reflects the oscilla¬ 
tion of gravitational waves. Our simulations have been terminated when the 
proper wave length Ly exceeds about 1.5/H in the contracting direction. It 
should be noted that Isaacson’s formula is not guaranteed if the proper wave 
length becomes comparable to the Hubble scale. Nevertheless, Fig. [3] shows 
that the maximum value of the deviation from Isaacson’s formula is about 
3% and the deviation of central value from the formula is even smaller. Next, 
we show the behavior of Lp^{H) and Li^{H). The behavior of LAo(Tf)(So in 
Fig. d]), Lai{H) (Si in Fig. [T]) and La 2 {H){S 2 in Fig. [T]) are almost same as 
each other. Therefore, we plot only Lao{H). 



Figure 4: The behavior of Lao (Ft)- Isaacson’s formula is ^AoiH) = yJbAo/H with 6 ao = 
0.224. 


11 







































As is shown in Fig. 01 Lxo{H) also oscillates in a similar way as Ly{H). The 
maximum value of the deviation from Isaacson’s formula around H ~ Hi q 
is about 4%, which is slightly larger than that of MH). 

The relation Li^{H) depends on the line position and can be classihed 
into six types (Fig. [5]). The behavior of and is similar to the 

behavior of Li^q{H) and 6i^o{H). Furthermore, the behavior of and 

Si^s^H) are almost same as and respectively, and the behaviors 

^L 2 (-ff) and have opposite phase to each other. Fig. [5] shows the value 

max(|(5Lo| : H ~ i^i.o) = 15% and max(|( 5 L 4 | : H ~ i^i.o) = 3%. On the other 
hand, max(|( 5 L 2 | : H ~ Ffi.o) = 40% and max(|(5L6| : H ~ Hi q) = 60%, which 
are the maximum in any 5i^{H ~ Hi q). 

In order to understand these behaviors and make the position dependence 
of proper length clear, we consider our model with the linear approximation. 
As is mentioned in the Sec. [ 2 ], our model corresponds to the linear standing 
gravitational waves. Since all lines are parallel to 2 ;-axis, only 7 ^^ can con¬ 
tribute to the values of proper lengths among components of the conformal 
metric. For the standing gravitational waves, the time evolution of '')zz is 
given by following form (see from appendix [B]) 



(28) 


where is a constant determined by the condition = 0. According 

to Eq. (1^ . on the line 0, 4 and 8 with coordinates (x,?/) = (0,0), (^, ^), 
(^, ^), 7 ^^ does not oscillate in the linear approximation. On the other 
hand, 7 ^^ on the line 2 and 6 with coordinates ( 0 , ^), (^, 0 ) oscillate with 
amplitude 2A. Thus, in our model, the position dependence of the effective 
scale factor L comes from the difference of the amplitude of gravitational 
waves on each line, which is affected by the constructive interference of two 
oscillation modes hH and (see Eq. flTT)) h This result implies that, when 
we use a proper length as the effective scale factor, position dependence must 
be carefully treated. 

4.2 The dependence of the initial amplitude A 

As is mentioned in the previous section, for the long wavelength (L > 1/hf), 
gravitational waves are superposition of the decaying mode and the growing 
mode and their ratio depends on the phase of gravitational waves at the hori¬ 
zon crossing. Due to this phase dependence, the behavior of gravitational 
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Llo{H): Isaacson’s formula is L^£q{H) = yJbijijH with = 0.224. 




Lli{H): Isaacson’s formula is L£f{H) = ^bLi/H with b^i = 0.223. 




Ll 2 {H): Isaacson’s formula is = ^/bL 2 jH with 6^2 = 0.220. 
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Ll 3 {H)-. Isaacson’s formula is L^£^{H) = ^JbL 3 /H with 6^3 = 0.223. 



Ll 4 :{H)-. Isaacson’s formula is L£^{H) = y/biAjH with 61,4 = 0.224. 




Llg{H): Isaacson’s formula is L£q{H) = with 64,6 = 0.218. 


Figure 5: The characteristic behavior of Ll{H) and Sl{H). 
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waves depends on A after the horizon crossing. In order to see this depen¬ 
dence, we focus on the time evolution of the Fourier component of with 
the wave number k = {27rnx/XoA'^'^^y/X q): 

lxx{t, k) = ^ d^xjxxit, x) cos(/c • x). (29) 

For /cq = (0, 0, 27r/Ao), the initial value of 7a;a;(^, ^o) is given by 

lxx{^^ ko) = — + 0{A^). (30) 

On the flat FLRW metric background, we can analytically calculate the 
behavior of linear gravitational waves. The metric can be expressed as fol¬ 
lows: 

ds"^ = —dt^ -I- a{tY {6ij + hij{t, F)} dx'‘dx\ (31) 

where a{t) is the scale factor of the background universe and hij(t,x) is the 
transverse traceless tensor of gravitational waves in the linear approximation. 
The evolution equation for a Fourier mode hij{t, k) of hij{t, x) is given by (see 
Appendix [B]) 


d \ 

k) + k) + —k‘^hij{t, k) = 0. (32) 

Behavior of hij{t, k^) for the short wavelength k^ > aH can be approximated 
by the WKB form as follows: 


hij{t, ko) ~ -hy 14=0 cos 

Qj 



(33) 


where we have used the normalization a| 4 =o = 1. Let us assume that the time 
evolution of the background metric is given by that of radiation dominated 
universe, 

a = A2bt + 1, (34) 

where b = by and a = Ly/Xo. From Eq. fl33l) . the time evolution of the 
Fourier component is given by 


a{'t)‘^hxxit, ko) = a— cos 


bXc 


(^V2bt + 1 — 


(35) 
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where we have set a?hxx\t=o = lxx\t=o up to the leading order and 0 is the 
integration constant. 

We want to evaluate the deviation between '"ixxit, ko) and a{t)‘^hxx{t, ko). 
Since we have used the WKB and the linear approximation to derive the 
expression fl3^ . when one of these approximations is violated, a(t)'^hxxit, ko) 
deviates from 'jxxityko). In our setting, this deviation may happen around 
the horizon crossing. After the horizon crossing even in linear regime, the 
“decaying mode” (for expanding universe) is enhanced and the deviation 
of our numerical solution from the linear WKB solution is expected. This 
deviation due to the decaying mode depends on the value of A. 

We introduce the deviation S{t) as follows 

“ a{t)^hxx{t, ko) 

2 - A^) -■ 

The evolution of 'jxxit, ko)/a(t), a(t)hxx{t, ko) and 6 (t) are shown in Fig. |6]for 
A = 0.1, A = 0.11 and A = 0.09. Around the horizon crossing, the behavior 
of 6 (t) depends on the value of A. So, one may expect that the deviation 
from Isaacson’s formula depends on A. Although the WKB approximation is 
violated around the horizon crossing, the values of 6^ and 5 a do not depend so 
much on the value of A. On the other hand, we can see significant dependence 
of ^ on 5 l- In particular, on the line 2 and 6, which correspond to anti¬ 
node of the standing waves, this dependence is large (for example, for A = 
0.07,0.08,0.09,0.1,0.11, max(|5L2(i7 ~ l/^v)|) is 60%, 40%, 60%, 30%, and 
40%, respectively). 

4.3 Temporal average 

As is mentioned in Sec. 13.21 to compare our results with the original Isaac¬ 
son’s formula, it is necessary to consider not only spatial average but also 
temporal one. The necessity of the temporal average is also explicitly shown 
in Appendix [B] for Isaacson’s formula in the FLRW universe. However, be¬ 
cause of the lack of the temporal periodicity, the temporal averaging is fairly 
ambiguous in the present case. We consider the following temporal averaging: 

^ rri+Xo/2 

{Ly)tem{r]) = ^ / dr]Ly, (37) 

^0 J»7-Ao/2 
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Figure 6: Evolution of the Fourier component 'yxx{t, ko)/a{t) (the left row) and the de¬ 
viation ((5S1) S{t) (the right row). The integral constant (j) in hxx{t,ko) is determined by 
fitting the numerical data in the region 0.0 < t < 1.0. We define the horizon crossing 
time tn and the big bang singularity time tp as = H{tH) ^ and L^^(to) = 0, 

respectively. 
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Figure 7: Evolution of dviri) = different grid sizes. The red, green 

and blue lines are A = Ao/80, A = Ap/lOO and A = Ao/120, respectively. The parameter 
b in the Ly = y/2bt + 1 is 0.223, which is given by fitting the numerical data in the region 
0.8 < Tv/A q < 1. 


where rj = dt/a is the conformal time. We note that the comoving wave¬ 
length of the linear gravitational waves is Aq. For simplicity we use the 
scale factor of the radiation dominated universe a = \J2bt -|- 1 □ as the scale 
factor in the dehnition of r], where b is given by htting \/2bt + 1 to Ly/\o 
in the region L^/Xq = [0.8,1.0]. Thus, the conformal time is rewritten as 
the following form: rj = a/b = 1/y/bH. The deviation by = ^-^v)tem0?)-ny(??) 

is depicted as a function of rj in Fig. 0 The value of 6y does not converge 
within our numerical precision because the 6y depends on A. Nevertheless, 
from this hgure, it is suggested that the deviation 6y is at most ~ 0.1%. 


5 Summary and discussion 

In this paper, we have investigated validity of Isaacson’s formula by solving 
the Einstein equations with technique of numerical relativity. By solving 
the Hamiltonian constraint, we numerically constructed the initial data of 
the universe which contains the nonlinear standing gravitational waves in a 
cubic box with the periodic boundary. Then the time evolution was per¬ 
formed with the COSMOS code based on the BSSN formalism. In order 
to investigate the validity of Isaacson’s formula, we calculated the relation 

^ In this discussion, we always consider expanding universe, thus the scale factor a is 
an increasing function of time. 
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between the effective scale factor and the Hubble parameter and compared 
it with Isaacson’s formula. The effective scale factors are dehned from the 
proper volume, proper area and proper length, and the Hubble parameter is 
defined from trace of the extrinsic curvature of the time slice. The results 
are summarized in Table [T] The deviation from Isaacson’s formula 6y is at 



L ~ 1/H 

L ~ 0.5/ H 

L ~ D.l/H 

volume 

<5vmax ~ 3% 

^Vmax ~ 2% 

<5vmax ~ 0.5% 

area 

^Amax ~ 4% 

(^Amax ~ 2% 

^Amax ~ 0.5% 

length 

(^Lmax ~ 60% 

<5Lmax ~ 30% 

(^Lmax ~ 5% 


Table 1: The largest value of the deviations from Isaacson’s formula are listed for each 
value of the effective scale factors L determined by the proper volume, area and length. 

most about 3% in the range of our numerical calculation given by L < 1/if. 
Although Isaacson’s formula is not guaranteed for L ~ 1/if, this formula 
has still a few % accuracy. Behavior of Lao{H), Lai{H) and La 2 {H) are 
similar to that of the proper volume. While, for the effective scale factor 
dehned by the proper length, deviation from the Isaacson’s formula depends 
on the line position. This means that, when we use the proper length as the 
effective scale factor, the position dependence must be carefully treated. We 
also discussed the temporal average of the effective scale factor given by the 
proper volume, which gives a central value of the temporal oscillation of the 
effective scale factor. Our calculation of the temporal average suggests that 
the deviation of the central value from Isaccson’s formula is at most ~ 0.1% 
in the range of our calculation. 

One might expect that, when the gravitational wavelength is much longer 
than the Hubble scale, our model approaches to the Kasner solution. How¬ 
ever, in our calculation, we could not proceed time evolution beyond L ~ 
1.5/H. At this time, Hamiltonian constraint is violated, and we could not see 
any typical feature of the Kasner solution. Even using the hner resolution, 
the violation time does not change. This time tc of the constraint violation 
is near the big bang singularity time to = ^ (see Eq. flMIB : 

0.3 

|to — tc| ~ (in the case of ^ = 0.1). (38) 

H 

Therefore, this constraint violation would be originated from the big bang 
singularity. We need more refined numerical method to analyze the struc- 
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ture of the spacetime around the singularity [TH]. Introducing another time 
coordinate, such as the e-folding number, we may investigate more details of 
the behavior in the early stage of the universe. We leave it as a future work. 
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Appendix 

A Time evolution equations 

Here, we summarize time evolution part of the Einstein equation in vacuum. 



-l-d' ^{—DiDja + 2Di In ^Dja + 2Dj In ^Dia 

+a[Rij - Rjij — 2DiDj In $ + ADi In A/Dj In ^ 

3 

+ ‘^{DkDHn^ -2Di,\n^DHn'^)^,,]}, (42) 


where a is the lapse function, f3 is the shift vector, is the Lie derivative 
with respect to /3, and the other variable dehned in Secj2j 
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B Isaacson’s formula in Friedmann universe 


Isaacson investigated the effective energy momentum tensor of low amplitude 
and high frequency gravitational waves [5l[6] in general background. In this 
appendix, we derive the effective energy-momentum tensor of low amplitude 
and high frequency gravitational waves in the spatially flat FLRW universe. 
The full metric can be written as 


-f a^(t) {5ij -I- T)} dx^dx\ (43) 

where a(t) is the scale factor and is low amplitude and high frequency 


perturbation which satishes the transverse and traceless gauge conditions: 

= 0 , d^hij = 0 . 

In order to consider the low amplitude and high frequency gravitational 
waves, we introduce a small parameter e -C 1 and assume the following 
order of perturbation: hjj ~ ^ ~ A// ~ e, where A is typical amplitude 
of gravitational waves A is typical wavelength of gravitational waves, I is a 
typical scales of background Hubble scale in the present case. Then, we find 
dkhij ~ e° and didkhij ~ e"h 

We decompose the Ricci tensor by the order of e as follows: 




(44) 


where we assign R® = 0{e^) for the background part. Each terms are 



dj - (aa+2d?)5ijy 


(46) 

(46) 

(47) 

(48) 

(49) 

(50) 




( 52 ) 


(51) 
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( 53 ) 


- ——d h - d h ■ 

+ —(9”/?”'(9 h - h ■ 

\ 2ci^ ^ ^ 2ci^ ^ ^ j^m'^ni’} 


- k^ h ■ 

IL y Lyyii 


and dot denotes the time derivative. 

In our ordering of the perturbation, the leading order equation is given 
by 


02 


d 


—hij{t,x) + ?,H—hij{t,x) - —V hij{t,x) = 0. 


dt^ 


Using the WKB approximation, the solution is 


hij{t,x) = - 


kd^k 

(2nf 


(.4<f«j(()e<« + c.c.), 


(54) 


(56) 


where u^{t) is defined by 

u^{t) = (56) 

I is the size of the comoving box and is an integration constant which 
satishes = 0 and = 0. 

In the next leading terms of the Einstein equations, we hnd the time 
evolution equation of the background. Since we are interested in the back- 
reaction effect of on the dynamics of the background metric we 
extract homogeneous and isotropic part of contribution from R^y. That is, 
we consider the following field equation 

(57) 

with the effective energy-momentum tensor is dehned by 

( 68 ) 

where () denotes a spatial average which smoothes out the inhomogeneities 
and extracts the spatially homogeneous and isotropic part. The average scale 
is assumed to be larger than the wavelength of gravitational waves and we 
explicitly describe how to take this average below. 

First, to preserve the covariance of each variables under the spatial ro¬ 
tation, we impose the following property for the average for any physical 
quantity Q(/c, k')\ 

{Q{k, k)u-,ug{t, x)) = {Q{k, k')Rk 
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( 59 ) 


{Q{k, -k))^^S{k + i-“, 

r 


and 


{Q{k,k')ulug{t,x)) = 

= {Q{k,k))^^^5{k-t). 

Substituting Eqs. fl55l) into Eq. flHTll and fl53ll . we obtain 


(60) 


(fiS) = 


(d?> = -( 


l^d^k (3k‘^ 


(27r)^ I 4 




2a4 




i-k) 


'V 


(61) 


k^ 


Pd^k f kik 

W) 

•^tk "^jk ' -^ik •^jk 




3kikj 
' 2a2 




+ 


.(fc) .(fc)* .(fc)* .(fc) 

■^ifc '^jfc “^ik “^jk 


The Ricci scalar is given by 


{RA = ^ / + 4f 

Then, is given by 


T, 


(GW) _ 


00 


flvk,^. 

Af 

^73 

cc 



(62) 


(63) 


(64) 


By extracting the trace part, we obtain 


qcw) 




Pd^k 


>1' 


(65) 


The hrst and the second terms inside the integration of Eq. (|65|) are tem¬ 
porally oscillating, and we may eliminate this term by taking a temporal 
average and get following form: 

PdPk 


y(GW) ^ 




121. , " - p 

Qo? J (27r)3 


A 


( 66 ) 
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According to fl64p and fl66p , we obtain traceless energy momentum tensor 
like the radiation fluid. It means that the energy density of the effective 
energy momentum tensor is proportional to 1/a^. 
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